% Plot the power at M2 (12.421 hours) of the observatory data on to the
% predicted power map


S = dir('/home/mnair/projects/tides/*_2010.mat');
S1 = dir('/home/mnair/projects/tides/*_2010_New.mat');
load '/home/mnair/projects/tides/M2_Abs_Colormap' map1;
open('/home/mnair/projects/tides/M2_Bz_Abs_0_720_1440.fig');
hold on
% periods = [4 4.8 6 8 11.967236 12  12.421 12.6583 23.934472 24 25.891];
% period_string = ['S6';'S5';'S4';'S3';'K2';'S2';'M2';'N2';'K1';'S1';'O1'];

%OK I don't need periods < 10

periods_12 = [11.967236 12  12.421 12.6583 ];
periods_12_text = [11.9 12.01  12.43 12.67 ];
period_string_12 = ['K2';'S2';'M2';'N2'];

periods_24 = [23.934472 24 25.891];
periods_24_text = [23.8 24.01 25.9];
period_string_24 = ['K1';'S1';'O1'];


for nstn = 1:length (S);
    
  eval(['load /home/mnair/projects/tides/' S(nstn).name ]);

  L = 1./(F*3600) >= 12.419 & 1./(F*3600) <= 12.4240;  
  
  if any(L)
  m2_amp = mean(sqrt(P(L,1)/(en-st)));
  
  % Find color to plot
  
  cidx = floor(63*(m2_amp)/(10)+1);
    
    if cidx > 64,
        cidx = 64;
    elseif cidx < 1,
        cidx = 1;
    end;

    figure(1);
    plot(longitude, 90-latitude, '.', 'color',map1(cidx,:),'MarkerSize',50);
  end;
    
  figure(2);
    subplot(121);
  L = 1./(F*3600) >= 11.5 & 1./(F*3600) <= 13;  
  plot(1./(F(L)*3600), log10(sqrt(P(L)/(en-st))), 'b', 'LineWidth', 2);
  hold on;
  text(periods_12_text,zeros([1,length(periods_12)]),period_string_12);
  a = axis;
  
  for kk = 1 : length(periods_12),
      line([periods_12(kk) periods_12(kk)], [a(3) a(4)],'LineStyle','--', 'color','k');
  end;
    
  title(['Amplitude spectra at ' S(nstn).name(1:3) sprintf(' (lat=%5.2f,lon=%5.2f,ndays=%7.1f,median amp=%5.2f)'...
      ,latitude,longitude,(en-st)/(24*60), median(sqrt(P/(en-st))))]);
  xlabel('Hours');
  ylabel('log(nT)');
  hold off;
  
  subplot(122);
  L = 1./(F*3600) >= 23.5 & 1./(F*3600) <= 26.5;  
  plot(1./(F(L)*3600), log10(sqrt(P(L)/(en-st))), 'b', 'LineWidth', 2);
  hold on;
  text(periods_24_text,zeros([1,length(periods_24)]),period_string_24);
  a = axis;
  
  for kk = 1 : length(periods_24),
      line([periods_24(kk) periods_24(kk)], [a(3) a(4)],'LineStyle','--', 'color','k');
  end;
    
  xlabel('Hours');
  ylabel('log(nT)');
  eval(['print -dpng  /home/mnair/projects/tides/' S(nstn).name 'power_spec' ]);
  hold off;
  pause;
end;